function V = slvar(X, w, vmean)
%SLVAR Compute the variances or weighted variances
%
% [ Syntax ]
%   - V = slvar(X)
%   - V = slvar(X, w)
%   - V = slvar(X, w, vmean)
%   - V = slvar(X, w, 0)
%
% [ Arguments ]
%   - X:        the sample matrix
%   - w:        the weights of the samples
%   - vmean:    the pre-computed mean vector of samples
%   - shared:   whether all dimensions share the same variance
%
% [ Description ]
%   - V = slvar(X) computes the sample variances for all elements of 
%     the vectors in X. 
%
%     Suppose X is a d x n matrix, with each column representing a 
%     d-dimensional vector. Then V will be a d-dimensional vector, 
%     with V(i) being the variance corresponding to the i-th value in 
%     the vectors.
%       
%   - V = slvar(X, w) computes the sample variances with weighted samples.
%
%     As in slmean, w can be given in either of the following forms:
%       - a 1 x n row vector giving the sample weights. 
%       - a K x n matrix giving k groups of weights over all samples
%       - empty, simply meaning not weighting.
%
%   - V = slvar(X, w, vmean) computes the sample variances with the mean
%     vector(s) of the samples given. In this case, the function will take
%     vmean as the mean vector without re-computing it.
%
%     If the samples are not weighted, or these is one group of weights,
%     then vmean should be a d x 1 vector, otherwise, it should be a
%     d x k matrix.
%
%   - V = slvar(X, w, 0) computes the sample variances considering the
%     samples have zero mean. 
%
% [ Remarks ]
%   - This function gives maximum likelihood estimates of the variance,
%     instead of the unbiased estimates. It means that in the non-weighted
%     version, it computes \sum (x-m)^2/n instead of \sum (x-m)^2/(n-1).
%
% [ History ]
%   - Created by Dahua Lin, on Oct 18, 2007
%   - Modified by Dahua Lin, on Dec 19, 2007
%       - add the functionality of computing on multiple groups of weights
%       - restructure the computation
%

%% parse and verify input arguments

assert(isnumeric(X) && ndims(X) == 2, ...
    'sltoolbox:slvar:invalidarg', ...
    'X be a 2D numeric matrix.');
[d, n] = size(X);

if nargin < 2
    w =[];
    k = 1;
else
    if isempty(w)
        k = 1;
    else        
        assert(isnumeric(w) && ndims(w) == 2 && size(w,2) == n, ...
            'sltoolbox:slvar:invalidarg', ...
            'w should be a k x n numeric matrix if it is non-empty.');
        k = size(w, 1);
    end
end

if nargin < 3
    vmean = [];
else
    if ~isempty(vmean)        
        assert(isnumeric(vmean) && ...
            (isequal(vmean, 0) || isequal(size(vmean), [d, k])), ...
            'sltoolbox:slvar:invalidarg', ...
            'vmean should be either 0 or a d x k matrix if it is non-empty.');
    end
end


%% compute

% compute mean vector

if isempty(vmean)
    vmean = slmean(X, w);
end

% prepare transposed normalize weights

if isempty(w)
    w = ones(n, 1) *  (1/n);
else
    if k == 1
        w = w' * (1 / sum(w));
    else
        w = bsxfun(@times, w, 1 ./ sum(w, 2))';
    end
end

% compute variances
if slisallzeros(vmean)
    V = (X .* X) * w;
else
    V = (X .* X) * w - 2 * (X * w) .* vmean;
    V = bsxfun(@plus, V, bsxfun(@times, vmean .* vmean, sum(w, 1)));
end
        






